library(Seurat)Attaching SeuratObjectlibrary(ggplot2)
library(ggpubr)EPCAM+analysis
Set important paths
dataDirs <- list(
day60 = list(
cellRanger = "/srv/gstore4users/p29781/o29804_CellRangerCount_2022-11-04--11-20-44/day60oticdiff_Library_418886_1",
seuratReport = "/srv/gstore4users/p29781/o29804_ScSeurat_2022-11-08--10-01-51/day60oticdiff_Library_418886_1_SCReport",
labeledObject = "/srv/gstore4users/p29781/additional-analyses/2022-11-28-relabel-day60"
),
day30 = list(
cellRanger = "/srv/gstore4users/p27889/o29467_CellRangerCount_2022-09-16--09-26-33/day30oticvesicle_Library_412321",
seuratReport = "/srv/gstore4users/p27889/o29467_ScSeurat_2022-09-22--15-04-32/day30oticvesicle_Library_412321_SCReport",
labeledObject = "/srv/gstore4users/p27889/additional-analyses/2022-10-04-relabel-day30oticvesicle"
),
day08 = list(
cellRanger = "/srv/gstore4users/p27889/o28146_CellRangerCount_2022-07-11--09-59-01/placodeday8test1",
seuratReport = "/srv/gstore4users/p27889/o28146_ScSeurat_2022-05-19--11-14-04/placodeday8test1_SCReport",
labeledObject = "/srv/gstore4users/p27889/additional-analyses/2022-08-19-Seurat-relabeled"
)
)
dataDirs$day60
$day60$cellRanger
[1] "/srv/gstore4users/p29781/o29804_CellRangerCount_2022-11-04--11-20-44/day60oticdiff_Library_418886_1"
$day60$seuratReport
[1] "/srv/gstore4users/p29781/o29804_ScSeurat_2022-11-08--10-01-51/day60oticdiff_Library_418886_1_SCReport"
$day60$labeledObject
[1] "/srv/gstore4users/p29781/additional-analyses/2022-11-28-relabel-day60"
$day30
$day30$cellRanger
[1] "/srv/gstore4users/p27889/o29467_CellRangerCount_2022-09-16--09-26-33/day30oticvesicle_Library_412321"
$day30$seuratReport
[1] "/srv/gstore4users/p27889/o29467_ScSeurat_2022-09-22--15-04-32/day30oticvesicle_Library_412321_SCReport"
$day30$labeledObject
[1] "/srv/gstore4users/p27889/additional-analyses/2022-10-04-relabel-day30oticvesicle"
$day08
$day08$cellRanger
[1] "/srv/gstore4users/p27889/o28146_CellRangerCount_2022-07-11--09-59-01/placodeday8test1"
$day08$seuratReport
[1] "/srv/gstore4users/p27889/o28146_ScSeurat_2022-05-19--11-14-04/placodeday8test1_SCReport"
$day08$labeledObject
[1] "/srv/gstore4users/p27889/additional-analyses/2022-08-19-Seurat-relabeled"resultDir <- "/scratch/jcarreno/sta426_project/results"Load the annotated object
anno <- readRDS(file = paste(dataDirs$day60$labeledObject, "/scData.rds",sep = ""))Original UMAP
DimPlot(anno, reduction = "umap", group.by = "ident")
hc <- c("ATOH1", "ANXA4", "GFI1", "CCER2", "POU4F3", "ATOH1",
"MYO7A",
"MYO6",
"PCP4",
"ESPN",
"TMPRSS3",
"BDNF",
"GNG8",
"POU4F3",
"OTOF",
"FSIP1",
"ZBBX",
"SKOR1",
"DNAH5",
"SCL26A5",
"GATA3",
"LMOD3",
"FGF8",
"TMPRSS3",
"INSM1",
"DNM3"
)Select EPCAM+ cells
anno.epcam <- subset(anno, cells = WhichCells(anno, expression = EPCAM > 1.5))
DimPlot(anno.epcam, reduction = "umap", group.by = "ident")
Recluster EPCAM+ cells with Louvain method
ElbowPlot(anno.epcam)
anno.epcam <- FindNeighbors(anno.epcam, dims = 1:19, k.param = 5)Computing nearest neighbor graphComputing SNNanno.epcam <- FindClusters(anno.epcam, algorithm = 1)Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 481
Number of edges: 3571
Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.8568
Number of communities: 14
Elapsed time: 0 secondsanno.epcam$sub_cluster <- as.character(Idents(anno.epcam))DimPlot(anno.epcam, reduction = "umap", label = TRUE)
DoHeatmap(anno.epcam, features = hc)Warning in DoHeatmap(anno.epcam, features = hc): The following features were
omitted as they were not found in the scale.data slot for the SCT assay: SCL26A5
Number of cells per cluster
cellCountperCluster <- data.frame(id = Idents(anno.epcam))
barplot = ggplot(data=cellCountperCluster, aes(x=id)) +
geom_bar() +
geom_text(stat='count', aes(label=..count..), vjust=-1) +
theme_minimal()
barplot + labs(x="Cluster", y = "Cell count")Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
ℹ Please use `after_stat(count)` instead.
UMAP for set of genes
Otic epithelium markers
oep <- c(
'EPCAM',
'CDH1',
'SOX2',
'SIX1',
'OC90',
'SOX10',
'FBXO2',
'LMX1A',
'PAX2',
'PAX8',
'DLX5',
'GBX2',
'JAG1',
'TBX2',
'COL9A2',
'OTOA',
'MYO5C',
'OTOL1',
'USH1C',
'PCDH9')for (marker in oep){
tryCatch({print(FeaturePlot(anno.epcam, reduction = "umap", features = marker))},
error=function(e){print("ERROR")},
warning=function(w){print("WARNING")})
}
#FeaturePlot(anno.epcam, reduction = "umap", features = oep)General Hair Cells Markers
ghc <- c('ATOH1',
'MYO7A',
'MYO6',
'PCP4',
'ANXA4',
'GFI1',
'ESPN',
'TMPRSS3',
'BDNF',
'CCER2',
'GNG8',
'POU4F3',
'OTOF',
'FSIP1',
'ZBBX',
'SKOR1',
'DNAH5',
'SCL26A5',
'GATA3',
'LMOD3',
'FGF8',
'TMPRSS3',
'INSM1',
'DNM3'
)for (marker in ghc){
tryCatch({print(FeaturePlot(anno.epcam, reduction = "umap", features = marker))},
error=function(e){print("ERROR")},
warning=function(w){print("WARNING")})
}
[1] "WARNING"
All other clusters
oc <- c('TTr',
'CLIC6',
'HTRC2',
'PLTP',
'CLDN3',
'TFAP2A',
'KRT8',
'TP63',
'KRT19',
'KRT5',
'CXCL14',
'DSP',
'KRT18',
'COL17A1',
'TYR',
'TYRP1',
'ECT',
'SOX10',
'EDNRB',
'POSTN',
'PPRRX1',
'TWIST1',
'VIM',
'PDGFRA',
'TWIST2',
'TTN',
'PAX7',
'ACTC1',
'MYLPF',
'TNNTI',
'TNNII')for (marker in oc){
tryCatch({print(FeaturePlot(anno.epcam, reduction = "umap", features = marker))},
error=function(e){print("ERROR")},
warning=function(w){print("WARNING")})
}[1] "WARNING"
[1] "WARNING"
[1] "WARNING"
[1] "WARNING"
[1] "WARNING"
[1] "WARNING"Heatmap for genes of interest
gi1 <- c('CDH1',
'SOX2',
'SIX1',
'OC90',
'SOX10',
'FBXO2',
'LMX1A',
'ATOH1',
'ANXA4',
'GFI1',
'CCER2',
'POU4F3',
'GATA3',
'MAFB',
'NEUROD1',
'TNTKR3',
'PRPH',
'NTRK3',
'STMN2',
'DCX',
'TUBB3',
'OTX1',
'ZIC1',
'ZIC2',
'PAX6',
'PAX3',
'OTX2')DoHeatmap(anno.epcam, features = gi1)Warning in DoHeatmap(anno.epcam, features = gi1): The following features were
omitted as they were not found in the scale.data slot for the SCT assay: TNTKR3
gi2 <- c(
'TTr',
'CLIC6',
'HTRC2',
'PLTP',
'CLDN3',
'TFAP2A',
'KRT8',
'TP63',
'KRT19',
'KRT5',
'CXCL14',
'DSP',
'KRT18',
'COL17A1',
'TYR',
'TYRP1',
'ECT',
'SOX10',
'EDNRB',
'POSTN',
'PPRRX1',
'TWIST1',
'VIM',
'PDGFRA',
'TWIST2',
'TTN',
'PAX7',
'ACTC1',
'TNNTI',
'TNNII'
)DoHeatmap(anno.epcam, features = gi2)Warning in DoHeatmap(anno.epcam, features = gi2): The following features were
omitted as they were not found in the scale.data slot for the SCT assay: TNNII,
TNNTI, PPRRX1, ECT, HTRC2, TTr
Expression dot plot
DotPlot(anno.epcam, features=gi1) + coord_flip()Warning in FetchData.Seurat(object = object, vars = features, cells = cells):
The following requested variables were not found: TNTKR3
DotPlot(anno.epcam, features=gi2) + coord_flip()Warning in FetchData.Seurat(object = object, vars = features, cells = cells):
The following requested variables were not found: TTr, HTRC2, ECT, PPRRX1,
TNNTI, TNNII
Show new cluster of HC in the original UMAP
cells.use <- WhichCells(anno.epcam, idents = '13')
anno <- SetIdent(anno, cells = cells.use, value = 'Reclustered HC')DimPlot(anno, reduction = "umap", group.by = "ident")
DimPlot(anno, reduction = "umap", group.by = "ident", cells.highlight = cells.use) + NoLegend()
Violin plots
All clusters
for (marker in oep){
print(VlnPlot(anno, marker))
}
for (marker in ghc){
tryCatch({print(VlnPlot(anno, marker))},
error=function(e){print("ERROR")},
warning=function(w){print("WARNING")})
}
[1] "ERROR"
Only EPCAM+ cells
for (marker in oep){
print(VlnPlot(anno.epcam, marker))
}
for (marker in ghc){
tryCatch({print(VlnPlot(anno.epcam, marker))},
error=function(e){print("ERROR")},
warning=function(w){print("WARNING")})
}
[1] "ERROR"
